function [Q_need, detail] = compute_svg_need(mpc, CB_steps, pf_lo, pf_hi, P_deadband)
%% 在“固定最佳 CB”的基础上，测算达到 PF≥pf_lo 所需的 SVG 给定 (Mvar)
% - pf_set = pf_lo（最省 SVG 的达标点）
% - SVG 允许双向（Q 可正可负）
% - 自估灵敏度 k = dQ_slack/dQ_svg，按 ΔQ 试探计算
% - 二分区间：[-Qbig,+Qbig]，Qbig ≈ |Qtarget-Qcurrent|/max(k,0.2) * 1.3

if nargin < 5, P_deadband = 0.00; end

%% 1) 固定“最佳 CB”（严格：整数台数 × 步长）
m = mpc;
Qstep_Mvar = m.cb(:,2) * m.baseMVA;     % 每步 Mvar（cb(:,2) 为 pu）
req_steps  = round(CB_steps(:));
Q_cb_applied = zeros(3,1);

keep_ratio = 0.05;   % 与阶段1一致
for i = 1:3
    bi    = m.cb(i,1);
    Qd0   = mpc.bus(bi,4);
    stepM = Qstep_Mvar(i);

    Qd_floor    = keep_ratio * max(Qd0, 0);
    Q_allow     = max(Qd0 - Qd_floor, 0);
    allow_steps = floor(Q_allow/stepM + 1e-9);
    appl_steps  = min(req_steps(i), allow_steps);

    Qcb = appl_steps * stepM;
    Q_cb_applied(i) = Qcb;
    m.bus(bi,4) = Qd0 - Qcb;
end

% 初始禁止 SVG 出力（稍后固定给定值）
assert(isfield(m,'svg') && isfield(m.svg,'gen_idx') && ~isempty(m.svg.gen_idx), ...
    'compute_svg_need: 未找到 mpc.svg.gen_idx。');
gidx = m.svg.gen_idx;
svc_bus = m.svg.bus;
m.gen(gidx,3) = 0;  m.gen(gidx,4) = 0;  m.gen(gidx,5) = 0;

% CB-only 潮流
res_cb = runpf(m);
assert(res_cb.success, 'compute_svg_need: CB-only 潮流不收敛');

% PCC
[Pcb, Qcb_pcc, ~] = slackPQPF(res_cb);
Pabs = abs(Pcb);

%% 2) 目标 Qtarget（最小达标：pf_set = pf_lo）
if Pabs <= P_deadband
    % 近零功率带：不强逼 PF，维持当前 Q
    Qtarget = Qcb_pcc;
else
    phi_set  = acos(max(min(pf_lo,1),0));        % 直接取下限（最省SVG）
    Q_set_abs= Pabs * tan(phi_set);

    % 沿当前 PCC 的无功方向（更省 SVG）；若 Q≈0，就选更近方向
    sgn = sign(Qcb_pcc);
    if sgn == 0
        cand = [+Q_set_abs, -Q_set_abs];
        [~,k] = min(abs(cand - Qcb_pcc));
        Qtarget = cand(k);
    else
        Qtarget = sgn * Q_set_abs;
    end
end

%% 3) 估计灵敏度 k ≈ dQ_slack/dQ_svg（用小扰动）
% 目的：把“PCC 目标缺额”换算成“SVG 给定需要的大致幅度”，用于设定二分区间
function [fval, res] = eval_f(Qsvc)
    mm = m;
    mm.gen(gidx,3) = Qsvc; mm.gen(gidx,4) = Qsvc; mm.gen(gidx,5) = Qsvc;
    Vm_here = res_cb.bus(res_cb.bus(:,1)==svc_bus, 8);
    if isempty(Vm_here), Vm_here = 1.00; end
    mm.gen(gidx,6) = Vm_here;
    res = runpf(mm);
    if ~res.success
        fval = NaN; return;
    end
    [~, Qslack] = slackPQPF(res);
    fval = Qslack - Qtarget;   % 目标方程
end

% 当前点的“目标误差”
[f0, ~] = eval_f(0);

% 用 ±ΔQ 试探灵敏度（取很小的 ΔQ）
dQ = max(0.05, 0.02*max(abs(Qcb_pcc),1));  % Mvar
[fL, ~] = eval_f(-dQ);
[fU, ~] = eval_f(+dQ);

% k_est ≈ ΔQ_slack/ΔQ_svg（用中心差分）
k_est = NaN;
if isfinite(fL) && isfinite(fU)
    % f(Q) = Q_slack(Q) - Qtarget ⇒ Q_slack 变化 = -(fU - fL)/2
    k_est = -(fU - fL) / (2*dQ);
end
% 容错：给个合理下限，避免除数过小
if ~isfinite(k_est) || k_est <= 0.2, k_est = 0.2; end
if k_est > 1.2, k_est = 1.2; end  % 一般不会超过 1

% 需要的 SVG 大小（粗估）
Q_deficit = -f0;                       % = Qtarget - Q_slack(Q=0)
Qsvg_est  = Q_deficit / k_est;

% 二分区间大小（×1.3 余量）
Qspan = 1.3 * max(0.5, abs(Qsvg_est));
Qbig  = max(0.5, Qspan);

%% 4) 区间 [-Qbig, +Qbig] 上二分
[fL2, resL] = eval_f(-Qbig);
[fU2, resU] = eval_f(+Qbig);

shrink = 0;
while (~isfinite(fL2) || ~isfinite(fU2)) && shrink < 6
    Qbig = Qbig/2;
    [fL2, resL] = eval_f(-Qbig);
    [fU2, resU] = eval_f(+Qbig);
    shrink = shrink + 1;
end

if ~isfinite(fL2) && ~isfinite(fU2)
    % 端点均不收敛：回退到CB-only
    Q_need  = 0; res_best = res_cb; hit = false;
    detail = pack_detail(); return;
end

% 若 0 点已满足（非常靠近目标）
if isfinite(f0) && abs(f0) <= 1e-3
    Q_need = 0; res_best = res_cb; hit = true;
    detail = pack_detail(); return;
end

% 是否跨零？若没有，取离零最近端点（容量/可行性限制）
if ~isfinite(fL2) || ~isfinite(fU2) || (fL2 * fU2 > 0)
    if ~isfinite(fL2) || (isfinite(fU2) && abs(fU2) <= abs(fL2))
        Q_need = +Qbig; res_best = resU; hit = false;
    else
        Q_need = -Qbig; res_best = resL; hit = false;
    end
else
    % 标准二分
    lo = -Qbig; hi = +Qbig; fl = fL2; res_best = resL;
    for it = 1:20
        mid = 0.5*(lo+hi);
        [fm, res_m] = eval_f(mid);
        if ~isfinite(fm), hi = mid; continue; end
        if abs(fm) <= 1e-3 || (hi-lo) <= 1e-3
            Q_need = mid; res_best = res_m; hit = true; break;
        end
        if fl*fm <= 0
            hi = mid;
        else
            lo = mid; fl = fm; res_best = res_m;
        end
    end
    if ~exist('Q_need','var'), Q_need = 0.5*(lo+hi); end
end

%% 5) 输出 detail
[P2, Q2, pf2] = slackPQPF(res_best);
detail = struct('Q_cb',Q_cb_applied, 'Qtarget',Qtarget, 'Qbig',Qbig, ...
                'k_est',k_est,'Qsvg_est',Qsvg_est, ...
                'Pcb',Pcb,'Qcb_pcc',Qcb_pcc,'P2',P2,'Q2',Q2,'pf2',pf2, ...
                'hit',logical(hit));
end

%% ====== 工具函数 ======
function [P_slack, Q_slack, pf_val] = slackPQPF(res)
slack_idx = (res.gen(:,1) == 1);
P_slack = sum(res.gen(slack_idx, 2));
Q_slack = sum(res.gen(slack_idx, 3));
S = hypot(P_slack, Q_slack);
if S < 1e-9, pf_val = 1.0; else, pf_val = P_slack / S; end
end
